##########################################################################################

library(data.table)
library(optparse)

##########################################################################################

option_list <- list(
    make_option(c("--data_type"), type = "character"),
    make_option(c("--clust"), type = "character"),
    make_option(c("--cluster_file"), type = "character"),
    make_option(c("--cor_t"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    data_type <- "exp"
    clust <- 12
    cluster_file <- "~/20231121_singleMuti/results/celltype_plot/mfuzz/pheatmap_exp.cluster.5.gene.tsv"
    out_path <- "~/20231121_singleMuti/results/celltype_plot/mfuzz"
    cor_t <- 0.5
}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

data_type <- opt$data_type
cluster_file <- opt$cluster_file
cor_t <- as.numeric(opt$cor_t)
out_path <- opt$out_path
clust <- as.numeric(opt$clust)

dir.create(out_path , recursive = T)

##########################################################################################

dat_mfuzz <- data.frame(fread(cluster_file))

set.seed(1234)

##########################################################################################
## 阳性转录因子
postive_tf <- subset(dat_mfuzz , cor > cor_t)$GeneExpressionMatrix_matchName
other_tf <- subset(dat_mfuzz , cor < cor_t)$GeneExpressionMatrix_matchName

## 计算富集
dat_mfuzz$TF <- ifelse( dat_mfuzz$gene %in% dat_mfuzz$GeneExpressionMatrix_matchName , "TRUE" , "FALSE" )
dat_mfuzz$postive_tf <- ""
dat_mfuzz$postive_tf <- ifelse( dat_mfuzz$gene %in% postive_tf , "positive" , dat_mfuzz$postive_tf )
dat_mfuzz$postive_tf <- ifelse( dat_mfuzz$gene %in% other_tf , "other" , dat_mfuzz$postive_tf )

tmp_dat <- table(dat_mfuzz$postive_tf , dat_mfuzz$group_name)[c("positive" , "other"),]
## 分每个cluster计算富集
result <- c()
for(clusN in colnames(tmp_dat)){
    a <- tmp_dat[,clusN]
    b <- apply(tmp_dat[,colnames(tmp_dat)!=clusN] , 1 , sum)
    tmp_fisher <- fisher.test(data.frame(a,b))

    p <- tmp_fisher$p.value
    er <- tmp_fisher$estimate
    L95 <- tmp_fisher$conf.int[1] 
    U95 <- tmp_fisher$conf.int[2]

    tmp_res <- data.frame( clusN = clusN , P = p , ER = er , L95 = L95 , U95 = U95  )
    result <- rbind(result , tmp_res)
}

out_name <- paste0(out_path , "/pheatmap_" , data_type , "cluster." , clust , ".motif_enrich." , cor_t , ".tsv")
write.table( result , out_name , row.names = F , quote = F , sep = "\t" )

out_name <- paste0(out_path , "/pheatmap_" , data_type , "cluster." , clust , ".motif_clust." , cor_t , ".tsv")
write.table( tmp_dat , out_name , row.names = T , quote = F , sep = "\t" )

